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Abstract 

An axisymmetric Meshless Local Petrov-Galerkin 
(MLPG) algorithm is presented for the potential and 
elasticity problems. In this algorithm the trial and test 
functions are chosen from different spaces. By a 
judicious choice of these functions, the integrals 
involved in the weak form can be restricted to a local 
neighborhood. This makes the method truly meshless. 
The MLPG algorithm is used to study various potential 
and elasticity problems for which exact solutions are 
available. The sensitivity and effectiveness of the 
MLPG algorithm to various parameters such as the 
weight functions, basis functions and support domain 
radius, etc. was studied. The MLPG algorithm yielded 
accurate solutions for all weight functions, basis 
functions and support domain radii considered for all of 
the problems studied. 


“mesh” or “cells” either to interpolate the solution 
variables or to evaluate the integrals that appear in the 
weak form, but rather needs only nodes that are 
randomly distributed in the domain. 

In the present work, MLPG algorithm is formulated and 
implemented for axisymmetric potential and elasticity 
problems. Various choices of the basis functions, forms 
of weight functions, and sizes of support domain are 
considered. The approximate solutions obtained by the 
MLPG algorithm are compared with the exact solutions 
for evaluating their accuracies. 

Axisymmetric Potential Problems 

Consider a Poisson’s equation for an axisymmetric 
problem bounded by a toroidal domain with its cross 
section defined by £2 as shown in Figure 1 , 


Introduction 


In the past decade, several meshless methods for the 
solution of partial differential equations in science and 
engineering appeared in the literature [1-4]. Meshless 
methods retain all the advantages of the finite element 
method without most of the disadvantages such as, 
element locking and discontinuous derivatives of the 
secondary variables across the element boundaries. 
Atluri and Zhu [4] presented a new and innovative 
meshless approach that uses Petrov-Galerkin weight 
functions instead of the traditional Galerkin weighted 
residual method. This meshless method is ‘truly’ 
meshless. This method, referred to as the Meshless 
Local Petrov-Galerkin (MLPG) method, does not need 


d 2 u 1 9 u d 2 u 

— H 1 r- 

dr 2 r dr dz 2 


g(r,z ) 


in £2 (1) 


with boundary conditions 

u = u on r u and q ~ q on F q (2) 


where T =r u + T q and q = du / dn . In a steady-state 
heat transfer problem, that is a typical Poisson’s 
problem, u is a temperature function, q is heat flux 
through the boundary, and g is function of internal heat 
generation. 

The solution for Eq. (1) is sought in a weighted residual 
manner as 
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2njcV 2 u- g)-vr-d rdz =0 (3) 

where v is a weight function. The factor 2 n in Eq. (3) 
comes from integration over the tangential coordinate 
8, because of axisymmetry. 
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Axisvmmetric Elasticity Problems 


(7) 


The stresses G,y and body force f for an elastic 
continuum need to satisfy the equations of equilibrium 
condition 

<r V j + f,= 0 (4) 

where comma indicates the partial differentiation with 
respect to the variable that follows the comma. 

The boundary conditions on V bounding the domain Q. 
are 


«(R)=5>( r >a 

k 

where w* is a fictitious value of u at node k and (fa (R) 
are the shape functions of node k. (Figure 2 illustrates a 
least squares approximation to the u k values at each of 
the randomly distributed nodes for a 1-D case.) The 
shape function, (fa, can be written as 

m . 

(P k = lPj(R)-lMRV'B(R)] jk (8) 

M 

with 


u = u on T u and t = t on T q (5) 


where u and t are displacement and traction vectors, 
respectively. 


The weighted residual statements for Eq. (4) for the 
axisymmetric case can be written as 


2*rf(^ + ^s-+i(cr r -a g ) + f r )- 

or oz r 

v-r-dr-d z = 0 


2 arj(- 

G 


dr„ da, 1 

IL. _i L J - 

r 

v, r dr-d z =0 


3 + -^ + - T n+ft)- 

or dz r 


( 6 ) 


A(R)=iw k (R)-p T (R k )-p(R k ) ( 9 ) 

k=\ 

B(R) = [w,(R) • p(R,),.. • w N (R) • p(R„ )] ( 10) 

where N is the number of the nodes in the domain of 
influence (see, for example, Ref. 4), w k ( R ) is a weight 
function defined later, and p T (R ) are the basis functions 
and are chosen to be 

p T (R) = l 1, r, z } (Ha) 

for linear representation (m = 3) and 

p r (R) = {l, r, z,r\rz,z 2 } (lib) 


where a n Oq, g z , and T rz are stresses in cylindrical 
coordinates ( r,6 \ z ) and v r and v^ are weight functions 
for r and z directions, respectively. In this study, the 
weight functions v r and v z are chosen to be identical. 

Since the MLPG algorithm for the potential and the 
elasticity problems are very similar, the potential 
problem is utilized to demonstrate the formulation in 
the following two sections. 

Moving Least Square Approximation 

In the present meshless method, a trial function, «(R), 
with moving least squares (MLS) scheme is used to 
approximate (in cylindrical coordinates (r, z); 

R = i-r + j‘Z, the radius vector from the origin) the 
distribution of the potential function u over a number of 
randomly distributed nodes as shown in Fig. 1. The 
trial function, u( R), can be expressed in a series form as 
[1-4] 


for quadratic representation (m = 6). 

Note that in Eqs. (9) and (10), the matrices A (R) and B 
(R) are of size m x m and mxN , respectively, and w k 
(R) is the weight function associated with the node k . 
The weight function w k (R) is chosen such that it is non- 
zero over an influence sub-domain surrounding the 
node k and zero outside of the influence domain [1-4]. 
In this study, the domain of influence of node k is 
assumed to be a circle of radius /*. The radius l k must be 
large enough to contain at least m nodes in each 
direction to prevent the numerical singularity of the A- 
matrix in Eq. (9). The weight function is chosen such 
that it is smooth and equals unity at the center of the 
sub-domain, and equals zero at the boundary and 
outside of the sub-domain. In this paper, three spline 
functions with C 1 , C 2 and C 3 continuity are used as 
weight functions: 

ForC 1 : 


w*(R) 


_ f 1 3 A + ^A 


(12a) 

A> 1 - 
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forC 2 
w*(R) = 

for C 3 : 
W *( R ) = 


l-10 A 3 + 15 A 4 -6 ft 5 
0 


0< A <1 
A>1 . 
(12b) 


l-35p/+84/o/-70p t 6 + 20p t 7 °^A^1 
0 


(12c) 

where p k = d k / l k is the normalized distance and 
^ is the distance between center of the 

support sub-domain R k and the general point R. The 
functions are shown in Figure 3a for a 1-D domain. For 
a 2-D domain, the weight function can be illustrated as 
a ‘smooth tent’ as shown in Fig. 3b for a Cl continuity 
spline function. 


MLPG Formulation 


In the current implementation of MLPG formulation, 
the weighted residual form of Eq. (3) is used. The 
essential boundary conditions are included in the 
weighted residual statement using a penalty method as 

2 n | J (V 2 m - g) • v ■ r * dr ■ dz 

(13) 

-a - J (M-M)-v-r-drj = 0 


where a is a penalty parameter which is chosen as a 
large number (lxl0 16 is used in this study.) Using the 
divergence theorem, Eq. (13) can be recast (dropping 
the constant 2n) as 



du dv du dv 
dr dr dz dz 


- r - dr * dz- ^ g 'V - r dr * dz 


+ 




•r- v dT 


- a • 


j(u-u)vrdr = 0 

r„ 


(14) 


where n r and n z are the direction cosines of the normal 
of the domain boundary along the r- and z-directions, 
respectively. 


Unlike the Galerkin method where the trial and test 
functions are chosen from same space, the Petrov- 
Galerkin method uses trial and test function from 
different spaces. In the present work, the trial functions 
are chosen as in Eq. (7), and the test functions are 
chosen so that the evaluation of integrals involved in 
Eq. (14) for the k th node is confined to a local 
neighborhood R k . This judicious choice of the test 
function is what makes this method truly meshless. In 
the present implementation, the test function, v(R,R*) , 

is defined in a similar manner as the weight functions in 
Eq. (12) by replacing the l k with l 0 i.e., 

v(R,R t ) = w*(R) =/(A) 0 < < 1 (15) 

= 0 p k > 1 , 

where p k = ||r -R t |// o is the normalized distance. 

For an internal node, l 0 is selected to be equal or less 
than the shortest distance between the node and the 
domain boundary (for example, nodes p and q in Figure 
4). For a node on the domain boundary, the sub- 
domain is a segment of a circular region (for example, 
node s in Figure 4). 

Substitution of the chosen trial and test functions in the 
weak form of Eq. (14) leads to 

[*]•{«} = {/} (16) 


where 


*,=j 


_ ff ^(R) 9i/(R,R t ) <Vj(R) 9v(R,R t ) 
dr dr dz dz 


rdrdz 


+a- J* i (R)-KR,R 4 )T-rfr 

■ v(R,R k ) r dr 


(17) 


and 


f k = jq(R)-v(R,R k )-r- dr 

T sq 

+a ■ j* H(R) v(R,R t ) r £?r 

r , u 

- J g(R) v(R,R k )-r-dr-dz 


( 18 ) 
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where q(R) is the prescribed value of normal flux on 

r . In Eqs. (17) and (18), the region £2 S is the support 

sq 

domain of v(R,R*) and is assumed to be a circle of 

radius l Q (see Figures 1 and 4). The T JM and n, are 
segments of the boundary formed by the intersection of 
Q s and T (see Figure 4 nodes s and q). The region Q s 
can be made as small as possible by a judicious choice 
of the test function, v. Numerical integration is used to 
evaluate the integrals involved in Eqs. (17) and (18). A 
12-point Gaussian quadrature numerical integration is 
used in the present study. 

The axisymmetric elastic problem is formulated on 
similar lines and hence this formulation is not presented 
here. 

Examples 

To evaluate the current MLPG algorithm for potential 
and elasticity problems, several patch test problems 
with exact solutions are considered. The sensitivity of 
the MLPG solution to the basis functions defined in Eq. 
(11) and the three weight functions w k (R) selected in 
Eq. (12) is studied. Since all of the potential and 
elasticity problems in this study are linear problems, the 
patch test models are defined with an arbitrary constant 
a. The various normalized radiuses of domain of 
influence ( l k /a ) and normalized radiuses of support 
domains (IJa) are used to study their influence on the 
accuracy of the results. 

The results of the MLPG method are compared to exact 
solutions. Two error norms (||e w ||i and \\e M \\ 2 ) are used 
to evaluate the effectiveness of various parameters. 
These norms are defined as 


K l (19b) 

where M is the total number of randomly distributed 
internal points in the domain at which the numerical 
solution is evaluated and compared to the exact 
solution. Note that these internal points are independent 
points and are not associated with the nodal points used 
in the models. A value of M = 50 is used in this study. 


Potential Problems 

Three potential patch tests involving Laplace- 
(V 2 w=0), and Poisson- (V 2 w = constant ) equation 
problems are considered. 

Figure 5 shows two models of the Laplace and Poisson 
equation problems used for the patch tests. These 
examples are analyzed using the MLPG algorithm for 
various prescribed boundary conditions for u and q on 
the boundaries. 

A triangular domain for two mixed boundary value 
problems for Laplace equation is studied. The exact 
solution for this patch test problem is u = r 2 - 2z~ + 3 z * 
The radius of the domain of influence l k is set to equal 
to 4a. Two different types of choices for l a are 
considered. In the first choice, the distance from each 
internal node to the closest boundary is computed. 
These distances are used as the values of l 0 for that 
node. Then the smallest value of all of the l 0 values is 
used for all nodes on the boundaries. With this 
definition, the support domain radius is different for 
different nodes. The MLPG algorithm with different 
values of l 0 for each node recovered the exact solution 
for the Laplace problems. In the second choice, the 
smallest value of all the l 0 values from the first choice is 
used for all the nodes. Once again, the exact solution 
was recovered. In these studies, the effects of the 
radius of the domain of influence, l k > were also 
evaluated by varying the radius from 4 a to 10 a. The 
accuracy of results are insensitive to the radius of the 
domain of influence. 

The sensitivity of the MLPG solution is studied by 
further varying the support domain radius l a . The 
convergence of the error norm is presented in Figure 6. 
All of the nodes in the model were set to have the same 
value of l 0 . The results show that, for Laplace equation, 
all three weight functions gave excellent solutions for 
both the mixed boundary values problems and for a 
variety of the normalized support domain radius, IJa. 

The problem involving a Poisson’s equation over an L- 
shaped domain is considered next. The exact solution 

for the patch test problem is u - r 2 + 2z 2 + 3z - 2 . 
The convergence of the solution for various values of l 0 
and for different boundary conditions is presented in 
Figure 7. Once again, the exact solution is recovered by 
MLPG algorithm for all idealizations demonstrating 
that the present MLPG algorithm is very efficient and 
accurate. 


7=1 


(19a) 


Finally, a heat transfer problem involving radial heat 
flow in a hollow circular cylinder with prescribed 
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constant temperatures Ti and T 2 at inner and outer 
surfaces, respectively, is considered (see Figure 8). The 
exact solution of the problem is 


T = T -(T - 7V . The model used in the 

ln(r 2 /r,) 

MLPG algorithm for the problem with 25 randomly 
distributed nodes and the boundary conditions used are 
shown in Figure 8. 


The exact solution for the temperature is a natural 
logarithmic function and since the exact solution cannot 
be represented by polynomial, a study of use of various 
basis functions (polynomial function) is undertaken. 
Three polynomial basis functions, linear, quadratic and 
cubic are considered. Note that the higher order basis 
function requires longer computer times. Figure 9 
presents the error norm for various basis functions with 
a C 1 weight function. The support domain radius, l 0i is 
set equal to 0.35 a for the nodes on the boundary and 
equal to the distance to the closest boundary for the 
internal nodes. Very accurate solutions are obtained for 
all cases studied. Figure 9 shows that the higher order 
basis function yielded more accurate results than the 
lower order basis function. 


Elasticity Problems 

Several elasticity problems were studied to evaluate the 
MLPG algorithm. The elastic material properties, 
Young’s modulus E = 100 GPa and Poisson’s ratio 
V =0.3, were used in these problems. As in the case of 
the potential problems, the accuracy of the results are 
insensitive to the radius of the domain of influence, 
if the value of l k greater than 4 a. 

First the axisymmetric elastic MLPG method was 
evaluated on problems involving axial tension and 
radial compression of a circular cylinder with different 
boundary conditions. Since exact solutions for these 
simple patch test problems can be represented by the 
polynomial functions, the meshless method, as 
expected, recovered the solutions to machine accuracy 
for all weight and basis functions all of the models. 

Lame’s Cylinder 

A hollow cylinder subjected to constant internal (p \ ) 
and external (p 2 ) pressures (Lame’s Problem) is 
considered (see Figure 10). For the problem analyzed, 
internal pressure pi, external pressure p 2 , the inner 
radius r 7 = a and the outer radius r 2 = 2 a were used. The 
models used for the problem are shown in Figure 10. 
The exact solution of the problem is 

<r r = C, - C 2 /r 2 =C, + C 2 /r 2 ;<r. =2-v-C, 


5 


1 + v 


U=- 


C,(l-2v)r + — 


u g = u= 0 ; 


where 


r _ Pi' r \ -Pz'h . r . PjZll r 2 r 2 
~ 2 2 ’ ^2 „2 „2 1 2 


( 20 ) 


r, -r. 


Figure 1 1 shows the variation of the error norm for the 
three different weight functions with the normalized 
domain size using quadratic basis function. All of the 
nodes in the model were assumed to have the same 
support domain radius, l 0 . The error norm reduces 
substantially as the domain size increases for all the 
three weight functions considered. Figure 11 also 
shows that the error norm of (^-weight functions is 
lower than those with C 2 -and C 3 - weight functions. One 
of the reasons why (^-weight function produces less 
error than the C 3 -weight function is possibly due to the 
fact that C 1 function has a larger “effective radius” than 
the C 3 function. The effective radius is defined to be a 
value of Pi in Eq. (12) as the value of the corresponding 
weight function w k (R) is too small to operative in the 
integrals in Eq. (17). For an arbitrary small threshold 
value of the weight function, \v k (R) = 0.01, the effective 
radius of C 1 function is 0.94 while the C 3 function is 
about 0.86 (see Figure 3(a)). This suggests that for 
Lame’s problem the effective radius appears to be 
influencing the accuracy of the results. 

An additional study of the domain size effect is 
conducted using a model with 25 randomly distributed 
nodes (Figure 12). Two different types of choices for l 0 
are considered namely equal and unequal size domain. 
In the first choice “equal size domain”, the l 0 is kept the 
same for all of the nodes and equals to the smallest 
distance from an internal node to the closest boundary 
(see Figure 12). In the second choice, “unequal size 
domain” is used. In this choice, for each internal node 
the l 0 is kept as large as possible and set to be the 
distance of the node to the closest boundary. For each 
boundary node, the l 0 is kept to the distance from that 
node to its adjacent boundary nodes (see Figure 12). 
The results in Figure 1 1 suggest that a more accurate 
solution (smaller error norm) is obtained when a larger 
domain size is used. Thus, the unequal size domain is 
expected to yield a more accurate solution. This 
expectation is confirmed in Figure 12 with the 
randomly distributed node model. 

The effects of the basis function are also evaluated, 
since the exact solution for the displacement in Lome’s 
problem can not be represented precisely by any of the 
three - linear, quadratic and cubic - polynomial basis 
functions used. The error norm obtained with each of 
the choices of the basis functions is presented in Figure 
13. The cubic basis function yielded the most accurate 
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results while the linear basis function produced a 
marginally better result than the quadratic function. 


Rotating Disk 


The next problem studied is that of a disk with uniform 
thickness rotating with a speed of co rad/sec. as shown 
in Figure 14. In this problem, all the boundaries are 
stress free while the disk is subjected to an inertial force 
per unit volume that equals to p -a) 2 r , where p is the 

density of the disk material and co is the angular 
velocity. The exact solution of the problem when the 
Poisson’s ratio is equal to zero is given by [see pages 
335-337, Ref. 5] 


<*r 


= ~po) 2 {r i l + r 2 2 ~ 
= ^pco 2 (r l 2 +r 2 2 + 
= 0; u z =0 




pco 2 




( 21 ) 


where again E is Young’s modulus of the material. 


Three models are used to study the convergence of the 
solution to nodal refinement (Figure 14). In the model- 
1, seven nodal points are placed in the r direction. In 
the model-2, the distance between adjacent nodes is 
half of the distance between adjacent nodes in model- 1 
(d=A/2.) In the model-3, the distance between adjacent 
nodes is a quarter of the distance in model- 1 (d=A/4). 
The number of nodes in z direction is kept the same for 
all of three models. 


Figure 15 shows the effect of the nodal arrangement 
with constant domain size (equals to the distance 
between the nodes of model-3, i.e. l 0 = A/4). The 
model-3 has most nodal points and produced, as 
expected, the most accurate solution. However, the 
model-3 requires about 4 times the CPU time of model- 
1. The effect of the basis functions is shown in Figure 
16 using model-2 with C 1 weight function. The cubic 
basis function yielded a more accurate result than both 
of the linear and the quadratic basis functions. 
However, the error norm of the linear basis function 
accuracy is about 0.1%. 


Figure 17 shows the effect of the domain size using 
model-2 and a cubic basis function. The results, once 
again, show that the large support domain radius yields 
a small error norm. 


Concluding Remarks 

A meshless Local Petrov-Galerkin (MLPG) algorithm 
is presented for analyzing potential and elasticity 
problems in axisymmetric domains. In this method, the 
trial and test functions are chosen from different spaces. 
By a judicious choice of these functions, the integrals 
involved in the weak form can be restricted to a local 
neighborhood. This makes the method truly meshless. 
The MLPG algorithm is implemented and its efficiency 
is studied with respect to three basis functions, three 
different weight functions, and different sizes of local 
support domain. 

Potential patch test problems involving Laplace and 
Poisson equations are used to evaluate the efficiencies 
of the algorithm. To evaluate the effectiveness of the 
MLPG algorithm for elasticity problems, several simple 
patch test problems, a Lame’s cylinder problem and 
rotating disk problem, are studied. 

The C ! weight function yielded more accurate results 
than C 2 and C 3 functions for the patch test problems 
studied. Since the C 1 function requires lesser computing 
time than C 2 and C 3 functions, it should be used in most 
problems. The larger value of the support domain 
radius, / 0 , tended to decrease the error norm. The 
variation of the radius of the domain of influence, /*, 
does not show a significant effect on the accuracy of the 
result. For most problems studied, the cubic basis 
function yielded most accurate results. Increasing the 
order of the basis function substantially increases the 
computing time, but reduces the error norm. 

In summary, the current MLPG algorithm yielded 
accurate solutions for all the potential and elasticity 
patch tests and other problems studied. A simple C 1 
weight function, with a large support domain radius and 
a simple linear basis function are recommended for 
accurate results. 
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Figure 1 A typical axisymmetric potential problem 



Figure 2 Approximation for u using MLS scheme for a 1-D case 
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(a) Various spline weight functions (b) C 1 spline weight function 

used ( 1 -D representation) (2-D representation) 

Figure 3 Spline weight functions used 



Figure 4 Convenient definitions of lo in the 
MLPG method 
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(b) L-shaped domain 

Figure 5 Models used in the patch tests of 
Laplace and Poisson’s Equations 
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0.01 0.1 1 0.01 0.1 1 
Support domain radius, (IJa) Support domain radius, (Ida) 


(a) Mixed boundary value problem -1 (b) Mixed boundary value problem -2 

Figure 6 Laplace problem involving a triangular domain - Convergence of the solution for 
various values of l 0 . ( V 2 w = 0 ,u exacl = r 2 - 2 z + 3 z Quadratic basis; /*= 4a) 



0.01 0.1 1 0.01 0.1 1 


Support domain radius, (/</«) Support domain radius, {Ida) 

(a) Mixed boundary value problem -1 (b) Mixed boundary value problem -2 

Figure 7 Poisson’s problem involving an L-shaped domain - Convergence of the solution for 
various values of l 0 . ( V 2 n = 8 , u exac , = r 2 + 2z 2 + 3z - 2 ; Quadratic basis; /*= 4a). 
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(a) Circular cylinder 

Figure 8 Heat transfer problem in a 
hollow circular cylinder 



Figure 9 Heat transfer problem involving a hollow 
circular cylinder - Effect of basis functions. 

(T = T -{T -T 2 )- — > C 1 weight function) 



Figure 10 Hollow cylinder subjected to internal 
and external pressures (Lame’s Problem) 


Support domain radius, (lja) 

Figure 11 Lome’s cylinder - Effects of domain 
radius and the weight functions 
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Figure 12 
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Equal Size ( IJa ) 


Unequal Size ( IJa ) 


Lame’s cylinder - Effects of support domain radius ( IJa ) - 
(Cubic basis function, C 1 weight function) 


P i 







Linear 


Quadratic 

Basis Functions Used 


Cubic 


Figure 13 Lome’s cylinder - Effects of basis functions (C 1 weight 
function; lo/a = 0.75) 
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3a 




Model -3 



Figure 14 Rotating disk problem and the Figure 15 Rotating disk - Effect of nodal refinement - 

three models used (C 1 weight function; Cubic basis function; lja= 0.125) 



Linear Quadratic Cubic 

Basis functions used 

Figure 16 Rotating disk - Effect of basis 
functions - (C 1 weight function; l a la = 0.25; 
model-2) 



Support domain radius, (lo/a) 

Figure 17 Rotating disk - Effect of support 
domain radius - (C 1 weight function; 
Cubic basis function, model-2) 
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